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ABSTRACT 


Based on theories of laminate anisotropic elasticity and interlaminar frac- 
ture, the complete solution structure associated with a composite delamination 
is determined. Fracture mechanics parameters characterizing the interlaminar 
crack behavior are defined from asymptotic stress solutions for delaminations 
with different crack-tip deformation configurations. A numerical method employ- 
ing singular finite elements is developed to study delaminations in fiber com- 
posites with any arbitrary combinations of lamination, material, geometric, and 
crack variables. The special finite elements include the exact delamination stress 
singularity in its formulation. The method is shown to be computationally accurate 
and efficient, and operationally simple. To illustrate the basic nature of com- 
posite delamination, solutions are shown for edge-delaminated [0/— 0/-0/0] and 
[(±0)/(±0)/9O°/9O°] s graphite-epoxy systems under uniform axial extension. Three- 
dimensional crack-tip stress intensity factors, associated energy release rates, 
and delamination crack-closure are determined for each individual case. The basic 
mechanics and mechanisms of composite delamination are studied , and fundamental 
characteristics unique to recently proposed tests for interlaminar fracture tough- 
ness of fiber composite laminates are examined. Influences of lamination, geome- 
tric, and crack variables on the delamination behavior are investigated. 


iii 



1. INTRODUCTION 


la an associate paper [1], the fundamental nature of stress singularity 
and asymptotic solution fields associated with a delamination between 
dissimilar anisotropic fiber composites have been studied, Lekhnitskii f s 
complex-variable stress potentials [2] in conjunction with an eigenfunction 
expansion method have been used in the formulation and establishment of the 
general solution. The eigenvalues, especially the ones which characterize the 
strength of stress singularity, for delaminations with different local 
traction boundary conditions near the crack tip have been determined. Of 
particular interest are the asymptotic deformation and stress governing the 
composite delamination fracture. Since logarithmic terras are absent in the 
solutions for the delamination stress and deformation [1], the general 
structure of the complete solutions consists of only a power-type 
eigenfunction series of both singular and higher-order terms. The set of 
unknown constants in the eigenfunction series solution for a delamination 
problem can only be determined by solving the complete boundary-value problem 
with a full consideration of overall composite geometry, lamination and 
material variables, remote boundary conditions, and end loading conditions. 

In a finite-dimensional fiber composite laminate with simple lamination 
variables and crack geometry, for example, a symmetric angle-ply [0/-0/-0/0J 
composite containing edge delaminations along the 0 and -0 ply interface, the 
complete laminate elasticity solution can be determined in an explicit form by 
the use of different analytical methods such as the boundary-collocation 
method [3], However, for a composite laminate having more than four plies, 
the aforementioned collocation method is not applicable, and a more general 
and versatile method of approach to the problem is needed. The situation 
could become extremely complex for delaminations in a composite containing a 
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large number of plies with different fiber orientations and laminar 
thicknesses. The almost unlimited number of variables in a general composite 
delamination problem requires the development of an advanced analytical method 
which can not only deal with the local singular behavior of the delamination 
but also take into account various combinations of lamination, material, 
geometric, and crack variables. Numerical methods such as finite element 
methods are considered to be most attractive because of their versatility in 
handling mechanics problems with complex structural geometry and material 
properties. In this paper, the currently developed laminate elasticity 
solution is incorporated in the formulation of a set of conforming finite 
elements with singular derivative fields [4]. The special crack-tip elements 
are shown later to be particularly suitable for modeling the composite 
delamination problem. The use of laminate elasticity solutions derived from 
the associated paper [1] permits the inclusion of exact delamination stress 
singularities in finite element formulation, leading to extremely accurate and 
efficient numerical solutions for studying the fundamental behavior of 
composite delamination with complex lamination variables and geometric 
parameters • 

Specific objectives of this paper are to: (1) establish complete 

solution structures for different delamination configurations to serve as a 
basis of formulating an advanced numerical method, (2) construct special crack 
elements of various kinds to model the composite delamination, (3) properly 
define interlaminar fracture mechanics parameters, e.g., stress intensity 
factors and energy release rates, for general composite laminates, (4) examine 
the fundamental behavior and associated characteristics of the composite 
delamination, and (5) assess influences of lamination, geometric and crack 
variables on the delamination response and composite failure modes. 
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In the next section, the general structure of complete solutions for 
delamination stress and displacement fields is given. Fracture mechanics 
parameters in terms of interlaminar crack-tip stress intensity factors and 
strain energy release rates are defined for various deformation modes of 
composite delamination. Special six-node, quasi three-dimensional crack 
elements for the present composite delamination problem are introduced in 
Section 3. Formulations of the singular elements and adjacent nonsingular 
eight-node isoparametric elements are briefly outlined. Solution strategy and 
computational scheme for the delamination problem, especially in the case that 
crack-surface closure (i.e. the contact problem) occurs, are discussed. 
Computational methods for evaluating stress intensity factors and energy 
release rates by using the singular finite elements are given in detail. In 
Section 4, solution accuracy and convergence are studied to demonstrate the 
efficiency and effectiveness of the present approach. The fundamental 
behavior and unique characteristics of composite delamination are examined in 
Section 5. Two graphite-epoxy laminate systems, i.e., symmetric angle-ply 
[0/-0/-0/Q] and symmetric [ (±9)/(±0)/9O°/9O° ] s composites, with delaminations 
emanating from laminate edges are studied. Delamination crack-tip deformation 
and fracture mechanics parameters are determined for each case. Influences of 
fiber orientation, ply thickness, and crack size on delamination failure 
mechanics and mechanisms are investigated also. 
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2. COMPLETE SOLUTION AND FRACTURE MECHANICS PARAMETERS 


2.1 Complete Solutions for Stress and Displacement 

As discussed in the associated paper [1], the eigenvalues 6 n can be 
determined by solving the transcendental characteristic equation Eq. 20 in 
[1], Depending upon the local crack-tip boundary conditions and the interface 
continuity conditions, the <5 n have the following values: 

(i) delamination with open crack surfaces 


fn -V2 


6 

n 


= < 


n 


(n - V 2 ) ± i^r 


(single root), 
(triple roots), 
(single root), 


(n = 0,1,2,...“); 


(ii) delamination with closed crack surfaces in frictionless contact 


(iii) 


5 =< 
n 

|"n -V 2 

l n 

(double roots), 
(quadruple roots) , 

(n 0,1,2,.. .°°) 1 

delamination with 

closed crack surfaces in 

frictional contact 


' n - V 2 

(single root), 


6 = < 
n 

n - 3 

(single root), 



n 

(quadruple roots ) , 

(n = 0,1,2,...°°). 


Once the values of S n are determined, the relationship among C^'s can be found 
and the complete solutions for displacement and stress can be established in 
explicit forms. For example, for a delamination with homogeneous local 
boundary conditions in a composite laminate subjected to planar loading, the 
complete solutions for stress and displacement components have the following 
expressions : 
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«i - H'ln j, (Rel ^Vr V2 l + I ”' b (k«)n A lk Z k n ' V2 l> 

n=0 k=l 


+ l ( [ c 2n b kf +c 3 n b ta ) l ^ < A ik Z k^ + [c 4n b ta ) > «V!» 


k=l 


<<=5„ IteI Ji i ta A ik^' 1/2+lY l + c 6n 1 " 1 j, “ta ) + 


oi 


(i = 1,2, 4, 5, 6), (1) 


u j ” ^ ( c in ^ {Re[b kn )r (j+3)k Z k ^ 1 + Im ^ b (k+3)n r (j+3)k Z k 2 1 > / Cn+ V 2 ) 

J n=0 k=l 


+ I ( [<=2n b kf +c 3„ b k' ) l Re l r ( j+ 3)k Z k +1/ <" +1 >l + b 4n b to ,I "‘l r O + 3)k Z k +1/ <'‘ +1 )l> 


k=l 


+ <C 5n ,tel J, b la r (j + 3)k Z k +V2+lT/(n+l/ 2 +il,)1 

iC"~ X 


+ C6nlmt k=l r (j + 3) kzk 


n+ V 2 +1Y/ (n+ l/ 2 +iy) ] } ) + Uq j (j = 1,2,3), (2) 


/ \ 

where b^ n ' are known eigenvectors corresponding to the unknowns c an for each 

given 5 n> and a oi and u Q j are known quantities from particular solutions for 

each individual case. We note that the are found to be the same for 

all 5 f s because of the involvement of the term e^ n< ^(4> = n and — rc) in the 
n 

0(6^) matrix. The constants and are defined as 
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A lk = p k’ A 2k= 

A 4k~ "V A 5k = W A 6k = "V 


lk = r 2k = n k> 

r 3k = w k» r 4k = p k’ r 5k = q k’ r 6k = t k» 

(3) 


where P k , r^, p k , q k , t k are related to ply stiffness matrix and can be found 
in [2,3], The expressions for stresses and displacements for delamination 
with closed crack surfaces, i.e., nonvanishing local traction boundary 
conditions, can be easily determined to have forms similar to Eqs. 1 and 2 but 
with slight modifications owing to different numbers of algebraic multiplicity 
of the eigenvalues involved. 

2.2 De lamination Stress Intensity Factors and Energy Release Rates 

The stress and displacement fields for a delamination are shown in the 
preceding section to possess general form of Eqs. 1 and 2 with unknown 
constants c Qn to be determined. A proper analytical or numerical method with 
the aid of global laminate boundary conditions and remote loading conditions 
is required to determine the detailed solution for the complete boundary - value 
problem. Since the interlaminar fracture is controlled by local stress and 
deformation, the asymptotic solution is of primary importance and interest in 
understanding the near-field behavior and fracture phenomenon. The asymptotic 
solution is recognized to be singular in nature and governed by the 
delamination stress singularities, which have been obtained in detailed in 
[1]. As pointed out in the associated paper [1], the singular eigenvalues 
depend upon the local delamination configuration; thus, distinct structures of 
asymptotic solutions are obtained for different crack-tip deformations. 

In general, the asymptotic solution for a delamination stress field can 
be written in the following form: 
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a. = 


i 

l 


3 

l 


[f 


j=l k-1 


(j) 

ik 


£ 


+ f 


(j) 

i(k+3) 


*k J 


(i = 1,2,3,... ,6), 


(4) 


where Z^ and Z^ have their origin located at the delamination tip. 

The and are functions of lamina material constants, lamination 

variables, and loading conditions. The integer £ is the total number of 
eigenvalues which satisfy the constraint condition, 

-1 < Re[6j ] < 0 (5) 

For the convenience of further developments, Eq. 4 may be rewritten as 


a ± “ 1 sp\r,<|>; <S.), (6) 

j-1 J 

where sp^ is the jth component (corresponding to the eigenvalue 6 j ) of the 
asymptotic stress in the polar coordinate system. 

In view of the asymptotic solution structures given in Eq. 4 and those 
given in [1,5,6], it is possible to define, in the context of interlaminar 
fracture mechanics, the delamination stress intensity factors and strain 
energy release rates in a manner consistent with those for a homogeneous crack 
and for the refined model of an interface crack between dissimilar isotropic 
media [7,8]. For example, in the case of a closed delamination with crack 
surfaces in frictionless contact, the stress intensity factors can be 
introduced by considering the crack-tip interlaminar stresses 02 , Og and 
(i.e., o y , T xy , T yz , or a ^ ^ , x r ^ , x z ^ ) along the ply interface <f> = 0 as 
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K t = lim /- 2 trr a 9 (r, 0 ) 
r+0 

( 7 a) 

K tt = lim, / 2 rrr oAr, 0 ) 
11 r>0 + 6 

( 7 b) 

Kj-jj = lim + /2irr o^(r,0) 
r*0 

( 7 c) 


where r = 0~ and 0 + denote the positions infinitesmally behind and ahead of 
the crack tip, respectively, because the normal traction a 2 is finite ahead of 
a closed crack tip where shear tractions and a 4 are singular. 

The strain energy release rate, G, and its components Gjl (i = 1 , 2 , 3 ) may 
be evaluated by using Irwin's virtual crack extension concept [ 9 ] as 
G = Gj + G2 + G^ 

Sci 

= lim 23a f ^ a 2^ r [u~ m \<5a-r ,tt) ~ u^ m+1 ^(5a-r,-Tr) ] 

6a+0 a 0 2 2 

+ o 6 (r,0)[uj m) (<5a-r,TT) - u^ m+1) (6a-r ,- tt) ] 

+ a ^(r , 0 ) [u^ m ^ ( 5 a-r , tt) - u^ m+1 ^(6a-r,-i:)]}dr, (8) 

where Sa is the length of virtual crack extension. The interlaminar stresses, 
°2» °6» and a 4* in E< 1* 8 can be obtained from the asymptotic stress field such 
as Eq. 4 . The corresponding displacements are also those of the asymptotic 
field equations discussed in the previous section. In terms of the 

lsminstion stress intensity solutions , the G and G^ for a. closed craclc 
can be shown to have a simple expression as 


G g 2 + G3 - j (A 2 K I]; + A 3 K xii ) 


( 9 ) 


In Eq. 9 , Gj is identically zero because of the displacement continuity across 
the closed delamination surfaces, 

u< m) <r,*> - -0 


-c < r < 0, 


(10a) 
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and A£ and A 3 are evaluated from [ 5 ] with 



- u< m+1) = A^/ 2 j 5a-r | 

( 10 b) 

u<”> 

- u^ +1 > = A ^/2 6 a-r 

( 10 c) 


as r -► 0 + . 

For a delamination with extremely small crack closure, the simplified 
model by taking the limiting case of a partially closed crack discussed in 
[5,6] is used. The stress intensity solutions Kjj and K IIL are the same as 
those given in Eqs. 7(b) and 7(c), but the Kj is defined at r * 0 + as 

K t = lim, /2irr o 9 (r,0). (11) 

r +0 2 

And the corresponding strain energy release rates G and Gj^ , then, have the 
form, 

G = G 1 + G 2 + G 3 = £ (A t K r + A 2 K n + A 3 K in ). ( 12 ) 

In the case that a fully opened delamination is assumed and the 
eigenfunction expansion series is used for the asymptotic solution, the 
delamination stress intensity factors are introduced as [ 3 ] 

3 —5 

K = lim I / 2 ir r j s^ J) (r, 0;6 ), ( 13 a) 

r +0 j=l 2 3 

K II = lim + E /2w r J s^ ; (r,0;6 ), (13b) 

r -*-0 j=l J 

3 ”6 

K IH = lim + l ft* r 3 s< j) (r,0;6 ). 
r +0 j=l 3 4 J 


(13c) 
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We note that in this case exact integration in Eq. 8 can be carried out and 
the strain energy release rate, G, can be determined easily. However, 
individual components can not be separated explicitly in the integration 
because of the complex mathematical structure of the asymptotic solutions. 
Other methods such as the hybrid singular finite element analysis [10,11] are 
needed to determine the values of individual G^. 
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3. COMPUTATION OF FIELD SOLUTIONS AND FRACTURE MECHANICS PARAMETERS 

We have thus far obtained the complete solutions for stress and 
displacement and the fracture mechanics parameters from the asymptotic 
solutions for each individual case. The unknown constants c an as well as 
and are to be determined by solving the complete boundary-value problem, 
involving remote loading conditions and global geometric and lamination 
parameters. Now consider a finite-width composite laminate containing 
delaminations under mechanical loading as shown in Fig. 1 in [1], As 
mentioned in Section 1, several numerical methods could be used to determine 
the complete solution for the delamination problem in a finite dimensional 
laminate. Owing to the singular nature of the interlaminar crack, the complex 
structural geometry, and numerous lamination parameters involved, the special 
numerical method employing recently introduced conforming singular finite 
elements is an attractive approach to the current problem. In this approach, 
the exact delamination stress singularities can be included in the formulation 
of the special elements. Thus, the unknown constants associated with the 
asymptotic solution and the fracture mechanics parameters governing the 
delamination behavior can be evaluated conveniently with a high degree of 
accuracy and a fast rate of solution convergence. 

3.1 Singular Crack— Tip Elements and Surrounding Nonsingular Elements 

In this study, we generalize the formulation and concepts of the 
conforming singular elements originally introduced for homogeneous isotropic 
elasticity problems [4] to the present quasi three-dimensional, anisotropic 
composite delamination problem. Formulation of the singular elements is based 
on selection of shape functions and their derivatives containing the exact 
eigenvalues which meet the constraint condition Eq. 5 of stress singularity 
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derived from the aforementioned eigenfunction analysis for the composite 
delamination. Detailed discussion of the nature of the singular elements has 
been given in [4]. Only relationships relevant to the current development are 
given here. 

Consider a six-node triangular element with three degrees of freedom per 
node and stress singularity at node 1 as shown in Fig, 1(a), The element has 
general rectilinear anisotropic material properties with elastic compliance 
j • By a proper transformation, any point in the element defined in global 
Cartesian coordinates (x,y,z) can be referred to both local polar coordinates 
(r,<{>,z) and triangular coordinates (p,£,z) with the origin located at node 
1. Within the element, the displacement components u are related to the nodal 
displacements q by the interpolation (or shape) function N as 


u = U(x,y) + u with U = N q, 

~ ~ ~o ~ ~ 

where u q are known quantities resulted from applied loading, and 


(14a,b) 


T 

q 


II 

{ u 1? 

u 2 * u 3 } ’ 

T / 

% = {u ol’ 

c 

o 

ro 

(1) 

(1) 

(2) (2) 

„(2) 

f ,(6) 

2 ’ 

q 3 ’ 

q l ’ q 2 ’ 

q 3 

q l 



N = N (p, 
~s _s 

6 J ) » 



(15a) 

(15b) 

(15c) 


in which the superscript T denotes transpose of the associated column vector, and 
the number in the superscript parentheses in Eq. 15(b) refers to the 
associated nodal number of q. Explicit expressions of the shape 
function N , involving proper eigenvalues 6* determined for each individual 
delamination problem from [1] and local coordinates, are given in Appendix 
1, Equation 14(b) can be written in a more explicit form as follows: 

6 .+1 

u = {p 3 MCO + L(p,?)} q, 


( 16 ) 
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where L(p,£) is the nonsingular polynomial part of N (p,£;<$ ), independent of 

~ s J 

the number and location of nodes along the 2-3 edge of the element. Using the 
minimum potential energy theorem and following the same procedure as the 
conventional displacement-based finite element formulation, we can construct 
the element stiffness matrix k for the crack-tip element as 


k = // B T CB dA, 

~s * ~s „ _s ’ 

A 

s 

where C is the material stiffness matrix, and B has the form 

- 

y 1 

B - 3(p J M + L). 

"S M A. 


(17) 


(18) 


In Eq. 18, the 9 is a matrix differential operator. We remark that the shape 
function of the singular element is chosen such that the element conforms with 
a nonsingular quadratic element matched through the common element boundary 
(i.e., along edge 2-3), and with singular elements of the same formulation 
through boundaries 1-2 and 1-3. 

The surrounding nonsingular elements used in this study are quasi three- 
dimensional, eight-node isoparametric elements [Fig. 1(b)] with 24 degrees of 
freedom (three D.O.F.'s per node). Formulation of the element stiffness 
matrix for the nonsingular element has been given in detail in [12,13]. The 
element stiffness matrix k f for the adjacent element can be shown as 


k 

.r 


where 


// B T C B dA, 
* ~r „ „r ’ 
A 

r 


(19) 


B = 8 N . 
~r . _r 


( 20 ) 


The shape function N r for an eight-node, quasi three-dimensional isoparametric 
element has standard quadratic expressions which can be found in [12,13]. 
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3.2 Solution Procedure and Iteration Scheme for Crack Closure 

In the finite element discretization of the delaminated composite, the 

singular elements are placed in a ring form with the delamination tip as the 

common node (Fig. 2). The element stiffness matrices k and k are true 

— s _r 

stiffness matrices relating unknown nodal displacements to nodal force 
vectors. The standard procedure of the matrix-displacement method [14] can be 
used to assemble the global stiffness matrix K and loading vector Q leading to 
the relationship. 


K q = Q + Q . 
~ ~ _ „o 


( 21 ) 


where Q q denotes the additional nodal force resulted from the applied strain 
e Q , and the assemblage may be expressed symbolically by 


K« I k (i) + l k ( J\ 


i=l 


j = l 


Q = 


s 

l 

i=l 


n 


Qg l) + l 0£ j) 

j = l ^ 


(22a) 


(22b) 


in which n s and n r are the total numbers of singular and nonsingular elements, 
respectively. 

In actual numerical computation of a delamination problem, a finite 
length of crack closure is assumed first. The problem now becomes an elastic 
contact problem because a part of the local boundary conditions is not known 
and needs to be determined from the solution. Specifically, the following 
continuity conditions are required along the closed portion of the 
delamination, -c < r < 0: 
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(a 2 (r,Tr)] = 0, (23a) 

lu 2 (r,ir)] = 0, (23b) 

where the bracket [ ] denotes the jump of the associated quantity across the 
closed crack surfaces, e.g., u ^ m+ ^(r,-if) - u^^r.ir) = 0. Thus, we solve 
Eq. 21 in terms of unknown pressures along the closed crack surface and use 
Eq. 23(a) to determine the contact stress on this surface. To enforce the 
conditions Eqs. 23(a) and 23(b), the solution technique proposed by 
Francavilla and Zienkiewicz [15] for an elastic contact problem is employed. 
The numerical procedure involves an iteration scheme to determine crack 
closure length and contact stress along the delamination surface. If the 
solution is admissible, the contact stress so obtained must be in compression 
and the displacement field should have no overlapping or interpenetration 
outside the contact region. In the case that crack closure length is found to 
be extremely small, say, less than the order of 10 ^ ~ 10 ^ times crack 
length, the delamination is then assumed to be open, and the simplified model 
discussed in [1,5] and in the preceding sections in this paper is used. For 
the case of a delamination with finite-length crack closure, the detailed 
iterative algorithm for evaluating crack closure and contact stress is given 
in Appendix 2. 


3.3 Conputation of Delamination Stress Intensity Factors and Energy Release 
Rates 

As mentioned in the preceding sections, the stress intensity factors and 

energy release rates for a delaraination are evaluated from the asymptotic 

solution of interlaminar stresses a 0 , a. and a, (or o , x and t ) and the 

2*6 4 ' y xy yz ' 

displacements u^ along the plane of the crack. For a finite dimensional 


composite laminate containing delaminations, the asymptotic stress and 
displacement can be conveniently determined by the aforementioned singular 
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finite element method. Along the delamination plane £ = £ 0 , the near-field 
stress and displacement fields are approximated by using Eq. 16 and its 
derivatives as 


o - p R(? 0 ;6) q, 


6 + l yj.f w V , 

U - P M(E ) q + u , 
~ „ o „ ~o 


(24a) 

(24b) 


where R(5 q ; 6) is a matrix of derivatives of the shape function M and ply 
elastic constants, and the p and £ 0 are related to the global coordinates by a 
simple transformation given in Appendix 1. 

For a partially closed delamination or a delamination with a very small 
size of crack-tip closure for which the simplified model with inverse square- 
root singularity is used, the stress intensity factors and strain energy 
release rates can be evaluated easily from the singular finite element 
results. Taking o and u along the delamination crack plane 

= 0 (i.e., 4> o = 0) , we can write the asymptotic interlaminar stresses and 
displacements along the interface in simple expressions as 


a 


i 



l/o 

u = B. r 1 + u , 
j k oj 


(i=2, k=l; i=6, k=2, and 
i=4, k=3) , (25a) 


(j=l, k=2; j=2, k=l , and 
j=3, k=3) , (25b) 


where A and B are obtained from the corresponding components of R(£ ;S) and M(£ ) in 
Eqs. 24(a) and 24(b) by setting = 0, and 6 = - V2 • Thus, the delamination 
stress intensity factors can be easily determined by 

K ± = fit A ± (i=l,II,III), (26) 

The energy release rates can be determined in a manner similar to that for 
through Eq. 8 as 
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G = Gj + G 2 + G 3 


lira 1 
6a +0 26 a 


6a 3 

/ { X ( \ 

0 k=l 



)[B k (6a-r) 



lim 1 
Sa+O 26a 


6a 

/ 


0 


X A k B k ( 

k-1 * k 


6a -r\^ 
r* ' 


dr. 


(27a) 


The term u Q j in Eq. 25(b) is not included, because it does not result in any 
contributions to G^ and G after integration. Integration of the singular 
integral Eq. 27(a) can be carried out explicitly without difficulty. Then the 
strain energy release rates have the form 

3 

G = l G i m 7T (A 1 B 1 + A 2 B 2 + A 3 B 3 ) * (27b) 

We remark here that each term in Eq. 27(b) corresponds to the individual G^ 
components and that for a delamination with finite crack closure, the first 
term in Eq. 27(b) is identically zero, i.e., Gj = 0, because of the continuity 
of displacement across the crack surface. 
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4. SOLUTION ACCURACY AND CONVERGENCE 

To demonstrate the accuracy and efficiency of the present method of 
approach, a symmetric angle-ply [45 0 /-45°/-45 0 /45°] graphite-epoxy laminate 
containing edge delaminations along the 45°/-45° ply interace is considered. 
For simplicity and without loss of generality, the composite is assumed to be 
subjected to uniform axial strain z z = z Q along the z-axis and have a geometry 
of b/h = 8 and = h^ = h = 0.25 inch and delaminations of length a = 0.25 
inch emanating from the edges. Elastic ply properties of unidirectional high- 
modulus graphite-epoxy identical to those of Eq. 28 in [1] are used. Owing to 
the geometric and lamination symmetry, only a quarter of the cross sectional 
area needs to be considered. In finite element discretization of the 
continuum, twelve special crack-tip elements of identical size and shape (Fig. 
2) are used to model the near-field response of the composite delamination. 

The crack-tip elements are embedded in the mesh of eight-node, quasi three- 
dimensional isoparametric elements. Local and overall mesh arrangements for 
the finite element analysis of the composite delamination are shown in Fig. 

3. To study the accuracy and convergence of field solutions, the mesh near 
the crack tip is continuously refined by halving the lengths of equal sides of 
the singular elements (e.g., the sides OE and OF in Fig. 3). 

Using the computational method and the solution scheme discussed in the 
preceding sections, numerical results are obtained for the convergence 
study. Significantly global crack closure with c/a = 0.34 is found for the 
delamination in the [45°/-45 0 /-45 o /45° ] graphite-epoxy under the uniform axial 
extension, resulting in a negative and identically vanishing G^. [The 
detailed nature of delamination closure and related problems will be discussed 
in the next section.] In Tables 1 and 2, delamination stress intensity 
factors and energy release rates associated with the finite-element mesh 
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refinement are shown for systemmatic reduction in crack-tip element size 
(Ah/h) and increase in the total number of elements. It can be seen clearly 
that stable and converged solutions for all and are obtainable, as the 
crack-tip element size becomes smaller than O.lh and the associated total 
number of elements exceeds 100. As anticipated, the delamination stress 
intensity factors (Table 1) are sensitive to the size of the crack-tip 
elements owing to the localized nature of the singular domain. However, the 
strain energy release rates are relatively insensitive to the element size and 
mesh refinement (Table 2) because the are related to the global structural 
response of the delaminated composite. Furthermore, the crack-closure length 
is also found to be insensitive to the mesh refinement in this case. 

To assess the accuracy of the current solutions for the composite 
delamination with global crack closure is not trivial because no analytical 
and numerical solutions are available in the literature for comparison. The 
only analytical study, which may be used as a reference, deals with the 
delamination problem by assuming the crack surface being fully open and 
employing an eigenfunction expansion method with the aid of a boundary 
collocation technique [3]. In Fig. 4, stress intensity solutions determined 
by the present singular finite element method including the crack-closure 
consideration are presented as a function of delamination length in the [45°/ 
-45°/-45°/45°] graphite-epoxy. The dominant stress intensity factor K m 
determined for the partially closed composite delamination by the present 
approach is about 6 ~ 7% higher than the value for an open crack from the 
boundary collocation results. However, solutions for Kj obtained by using the 
two distinct models differ from each other both in sign and in magnitude. In 
the current singular finite element analysis of the delamination problem, Kj 
is found always to be negative owing to the aforementioned crack closure. 
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whereas Kj determined for an assumed open crack is positive. Values of Kjj 

and G 2 determined by the two approached are both negligibly small. The strain 

energy release rates obtained for the delamination are given in Fig. 5. For 

the partially closed delamination determined by the present approach, is 

62 

identically zero, whereas the value of G^/CIO e^) from the assumed open crack 
model is of the order of 10 -6 — a value much smaller than that of G 3 , the 
dominant component of energy release rates. Despite the different nature at 
the crack tip in these two approaches, the values of G 3 differ only slightly 
from each other. 

We remark that though mathematically rigorous, the solution obtained 
using the eigenfunction expansion method with the assumption of open crack 
surfaces is physically inadmissible because severe interpenetration of 
dissimilar materials is found in the oscillatory displacement solution near 
the crack tip. Thus, the present singular finite element approach with a 
partially closed crack consideration provides physically meaningful solutions 
with a high degree of efficiency and accuracy for general composite 


delamination problems. 
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5. THE FUNDAMENTAL BEHAVIOR AND CHARACTERISTICS OF COMPOSITE DELAMINATION 

In the preceding sections, we have introduced the complete solution 
structure and defined the governing fracture mechanics parameters for the 
delamination problem. Special finite elements with singular derivatives have 
been formulated to model the near-field response of a delamination in a finite 
dimensional composite laminate. Solution convergence and accuracy have been 
affirmed in terms of computational parameters (e.g., the size of crack-tip 
elements, the degree of finite element discretization, etc.). Having 
established these basic mechanics theories and the numerical method, we now 
proceed to study the fundamental behavior and characteristics of delamination 
in fiber composites with general lamination and geometric variables. 

Two fiber composite material systems are examined in this section: (1) 

symmetric angle-ply [±0] s graphite-epoxy laminates with edge delaminations 
between 0 and -0 plies, and (2) symmetric [(±0)/(±0)/9O°/9O°] s graphite-epoxy 
with edge delaminations between -0 and 90° plies. The symmetric angle-ply 
composite laminate system is selected because several unique delamination 
characteristics are observed, which can be used to illustrate most clearly the 
basic interlaminar fracture mechanics and failure modes. More importantly, 
some of the most fundamental nature of delamination fracture determined from 
the presently introduced physically admissible model and mathematically 
rigorous solutions are not observable in the previously obtained solutions 
which contain the inadmissible oscillatory stress and deformation [3,5]. The 
[ (±0)/ (±0)/9O°/W^] s graphite-epoxy system is studied also because this 
lamination system, especially the one with [ (±30° ) / (±30° )/90°/90°] s fiber 
orientations, is currently being considered for use in the evaluation of 
interlaminar fracture toughness of composite materials under static and cyclic 
loading [16,17]. The ply elastic properties of high-modulus unidirectional 
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graphite-epoxy given in Eq. 28 of Ref, [1] are used in all computations. In 
the [0/-0/-0/0] composites, the following lamination and geometric parameters 
are employed: hi=h2=h=0.25 inch and b/h = 8, whereas in the 

[(±0)/(±0)/9O°/9O°] s composite laminates, actual dimensions of unidirectional 
graphite-epoxy used in laboratory experiments [16,17] are taken, i.e., 
hi = h 2 ■ ... =» hjj = 0.0054 inch and b = 0.75 inch. Moreover, in all cases 
studied in this paper the composite laminates are considered to be subjected 
to uniform axial extension e z = e Q along the z direction. 

Since the interlaminar crack-tip deformation and fracture mechanisms are 
governed by the asymptotic field solutions , we shall examine the fundamental 
behavior and associated characteristics of composite delamination in terms of 
interlaminar fracture mechanics parameters, i.e., crack-tip stress 
intensities, strain energy release rates, and crack-surface closure. 

5.1 Influence of Fiber Orientation 

The behavior of a delamination in the relatively simple angle-ply [±0] g 
laminates is significantly influenced by the fiber orientation 0. Results 
obtained by using the iterative solution scheme in Section 3.2 reveal that the 
delamination always possesses a finite-length crack-tip closure. Assuming 
that the crack surface is in frictionless contact, we find that the closure 
length is global in general. For example, in the [±0] s graphite-epoxy with 
edge delaminations of length a = 0.25 inch, significant crack closure is found 
in each case studied (Table 3). The contact-zone size, c/a , varies from 
approximately two-tenths to more than one-third of total delamination length 
for 0 ranging between 15° to 60°. The crack-tip closure results In a negative 
opening-mode stress intensity factor Kj < 0 and an identically vanishing 
energy release rate Gj = 0, as shown in Tables 4 and 5. Thus, the 
delamination behavior in the [0/-0/-0/0] composites is apparently governed by 
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interlaminar shear stresses and deformations. In Tables 4 and 5, the out-of- 
plane tearing-mode (mode III) stress intensity factor K-j-j-j- and associated 
energy release rate G 3 are clearly the dominant ones controlling the 
delamination fracture, (The values of and G 3 are several orders of 

magnitude higher than those of Kj-j- and G 2 .) Also, the composite with 0 - 15° 
possesses the highest Kj-q and G 3 among all 9's studied. These results 
suggest that in symmetric angle-ply composites, delamination initiation and 
growth are more intimately associated with local interlaminar shear than 
transverse normal stress. 

We reiterate that the solutions obtained here by the use of the present 
interlaminar crack-closure model and the rigorous mathematical (combined theo- 
retical and numerical) approach are physically admissible and meaningful, as 
contrary to the previously obtained inadmissible, oscillatory solutions from 
an eigenfunction expansion approach with the assumption of a fully open crack. 

In the more complex graphite-epoxy laminates with [(±0)/(±0)/9O°/9O°] o 

s 

ply orientations, delaminations are always observed to occur between the -0 
and 90° plies [16,17], Starting with the partially closed crack model and 
using the iterative solution scheme, we find that the closure length of a 
delamination is less than 10 “ 6 inch. Thus, for this extremely small crack 
closure the simplified model with an inverse square-root stress singularity 
discussed in Section 4.3 of [1] and in [5,6] is employed. In Figs, 6 and 7, 
variations of stress intensities and energy release rates G ± with fiber 
orientation 0 are shown for a delamination of length a = 0.5b, The crack tip 
is apparently governed by the opening- and inplane shearing-mode stresses; the 
Kj and K-q are one order of magnitude higher than for all 0 's studied. 

The value of G 3 is vanishingly small in general, and Gj and 63 are three 
orders of magnitude higher than G 3 because of the extremely small value of 
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A-j. Note that values of G^, G2 and G3 depend not only on but also on as 
indicated in Eqs. 12 and 27(b). [Unfortunately, an explicit relationship 
among A^, ply elastic constants, lamination parameters, geometric, and loading 
variables needs to be determined numerically because of the complex algebra 
involved.] We remark that the maximum Kj and Gj appear at 9 = 32.5° and the 
maximum Kjj and G 2 occur at 9 = 35° in the [(±9)/(±9)/90°/90°] s graphite-epoxy 
system. The selection of [ (±30 o )/(±30°)/90 o /90° ] s lamination in the recently 
proposed interlaminar fracture toughness tests [16,17] using edge -delamination 
composite specimens is therefore a proper choice. 

5.2 Influence of Ply Thickness 

The ply thickness is an important geometric variable in studying the 
delamination behavior of composites. Changing the laminar thickness in a 
laminate alters the lateral constraint of adjacent plies in the thickness 
direction and, thus, directly affects the interlaminar crack behavior. For 
illustration, fracture mechanics solutions for the [45 0 /-45°/-45 0 /45°] 
graphite-epoxy composites with several ply thicknesses are presented in this 
section. For simplicity and without introducing further complications, the 
laminate width 2b and thickness 2W as well as the delamination length a are 
kept constant as before , while variation of ply thicknesses hj^ and h 2 (with hj 
+ h 2 = W = 0.5 inch) is considered. 

As shown in Figs. 8-10, altering the ply thickness t^/W has appreciable 
effects on the fundamental behavior of the delamination. For instance, when a 
thick outer 45° ply is used, say, hj/W > 0.8, the entire delamination surface 
is in contact with the other (Fig. 10), leading to a negative K-j- with Gj = 

0. Delamination fracture in this situation is governed by the tearing-mode 
stress intensity factor Kjj-j- and the associated energy release rate G3 (Figs. 

8 and 9). When equal ply thickness, i.e., hj = h 2 , is used, the in-plane 
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shearing-mode stress intensity and the energy release rate G 2 are about 
zero, and Kjjl has the maximum stress intensification. However, as the outer 
45° ply becomes thin, e.g., hj < 0.2W, the surfaces of delamination are opened 
up, i.e., c/a = 0, as shown in Fig. 10. In this situation, contributions of 
Kj and Kjj to the local failure become increasingly appreciable; Kjj and K m 
are of the same order of magnitude but G 2 and Gj are still relatively small as 
compared to Gg. In the limiting case of h^/W ■* 0, the in-plane shear could 
become dominant to govern the delamination. 

5.3 Influence of Interlaminar Crack Length 

For the symmetric composite laminates with either [±9] g or 
[(±0)/(±0)/9O°/9O°] s fiber orientations, the interlaminar crack-tip behavior 
is affected by the size of delamination. As discussed in Section 5.1, the 
delamination in the [0/-9/-0/0] graphite-epoxy system has a global crack 
closure. Figure 11 reveals general characteristics of the crack closure in 
the delaminated [45 0 /-45°/-45°/45 0 ] graphite-epoxy. Under a uniform axial 
strain e Q , the crack closure increases raonotonically with crack length until 
the delamination becomes about two-ply thicknesses. As the delamination 
extends further, the closure length approaches an asymptotic value of c/h * 
0.375, indicating that crack growth is governed by the interlaminar shear 

yz 

and the tearing-mode stress intensity factor Kjjj, as shown in Fig. 4. The 
values of T X y, Kjj and G 2 are orders of magnitude smaller than those of Xy 2 , 
Km anc * ^3 (Figs. 4 and 5). Note again that crack closure occurs for all 
a/h's in the [45 0 /-45°/-45 0 /45° ] graphite-epoxy, leading to an identically 
zero Gj and negative opening-mode stress intensity factor K-j- < 0. In the case 
of a very small delamination emanating from a free edge under rising load, the 
interlaminar crack is inherently unstable and extends rapidly to about one or 
two-ply thicknesses (i.e., the plateaus in Figs. 4 and 5) before stable growth 
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would occur. This phenomenon has been indeed observed in laboratory 
experiments of edge-delamination tests of [(±30°]/±(30 o )/90 o /90°] s T300-5208 
graphite-epoxy composite laminates [16]. 

In the [(ir30°)/(±30 o )/90 o /90°] s graphite-epoxy containing delaminations 
between -30° and 90° plies, the delamination is found to be open all along the 
crack surface for any given &/h, as discussed in Section 5.1. The solutions 
for and as a function of delamination length a / b are given in Figs. 12 
and 13. Also shown in Fig. 13 are the results reported in [16] by using a 
conventional nonsingular finite element method in conjunction with a virtual 
crack-closure scheme for calculating the strain energy rates [18]. The 
difference is very appreciable to warrant the necessity of using the advanced 
analytical and numerical technique for solving the delamination problem 
accurately. We observe that local delamination growth in the 
[(±30 o )/(±30°)/90 o /90°] s graphite-epoxy is governed by an inherently three- 
dimensional mixed-mode fracture process because the simultaneous presence of 
significantly high values of K-j- , K^x and Kxxi shown in Fig. 12. Owing to the 
vanishingly small A 3 , the value of G 3 is found to be negligible in comparison 
with Gj and G 2 (Fig. 13). We remark that in the limiting case of a very small 
delamination crack (i.e., a/h + 0 ), interactions occur between the 
delamination crack tip and the laminate edge, leading to a slightly higher 
value of &xil as shown in Fig. 12. 

A salient feature shown in Figs. 4, 5 and 12, 13 is that the stress 
intensity factors and energy release rates G^ become independent of the 
crack length a/h or a/h as the delamination extends beyond a few ply 
thicknesses. This unique feature has been observed and used in the 
experimental study of interlaminar fracture toughness of composite laminates 
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[16,17]. Provided that the and G.^ are accurately determined by appropriate 
methods such as the present one, the stress intensity factors and energy 
release rates associated with observed edge delamination initiation and growth 
may be useful for characterizing delamination fracture and interlaminar crack 
resistance of composite laminates. 
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6. SUMMARY AND CONCLUSIONS 

General expressions of stress and deformation fields in a fiber composite 
laminate with delaminations have been obtained from presently developed 
anisotropic laminate elasticity theory and basic concepts of interlaminar 
fracture mechanics. Fracture mechanics parameters such as stress intensity 
factors and energy release rates for delaminations with different crack-tip 
deformation configurations are defined in a manner consistent with those for a 
homogeneous crack and for the refined model of an interface crack between 
dissimilar isotropic solids [7,8], In finite-dimensional composite laminates 
with complex lamination and geometric variables, an advanced numerical method 
employing special singular crack-tip finite elements is developed for modeling 
the delamination. Exact delamination stress singularities obtained from the 
laminate elasticity solution are included in the crack-tip element 
formulation. Solution convergence and accuracy have been studied to ensure 
the validity of the results and to demonstrate the efficiency and 
effectiveness of the method. To illustrate the fundamental nature of 
composite delamination, numerical results are shown for the [9/-0/-0/8] and 
[ (±0)/(±0)/9O°/9O° ] s graphite-epoxy laminate systems containing edge 
delaminations under uniform axial extension. Fracture mechanics parameters 
and failure modes associated with the composite delamination are determined 
for each case. The basic mechanics and mechanisms of delamination are studied 
for the composites with different lamination and geometric variables and crack 
parameters. Based on the results obtained, the following conclusions may be 


reached : 
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(1) To study the fundamental mechanics of delamination in a composite 
laminate having a large number of plies with different fiber orientations, 
advanced numerical methods such as the present one are essential to take into 
account the crack-tip singularity and the large number of lamination, 
material, and geometric variables involved. 

(2) In the crack-tip element formulation, the inclusion of delamination 
stress singularities determined from the laminate elasticity solution by using 
an eigenfunction expansion method leads to very accurate solutions with a 
rapid rate of convergence. This is particularly advantageous for studying 
delamination problems in composite laminates with complex lamination, 
geometric, and material variables. 

(3) Since the singular eigenvalues and their algebraic multiplicity are 
different for different crack tip deformation configurations, the crack 
element formulation and solution strategy for each delamination problem need 
to be treated on an individual basis. Stress intensity solutions and energy 
release rates in each case should be evaluated in accordance with the 
appropriate delamination models given in Section 2. 

(4) The state of stress and deformation in the vicinity of a delamination 
crack tip are three dimensional in general. The asymptotic solutions can not 
be determined accurately by using classical lamination theory nor by any 
approximate methods without including interlaminar stresses and the correct 
stress singularities associated with the delamination. The current laminate 
elasticity solution and associated numerical method provide accurate 
information on the singular nature of the crack tip and complete field 
solutions for the delamination problem. 
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(5) The singular delamination stress field may be characterized by the 
presently introduced interlaminar fracture mechanics parameters. Crack-tip 
stress intensity factors and energy release rates G^ can be determined only 
by solving the complete boundary-value problem. In contrast to the singular 
eigenvalues which depend only on local boundary conditions and material 
properties of adjacent plies, the and G^ are functions of all lamination 
and geometric parameters, remote loading conditions, and crack variables. 

(6) Owing to the local nature of interlaminar fracture, the mechanics and 
mechanisms of delaraination growth are governed by the crack-tip stress 
intensity factors. In angle-ply [0/-0/-G/9] graphite-epoxy, crack-tip closure 
occurs and, thus, < 0; delamination growth may be, therefore, more 
intimately related to interlaminar shear stresses and interface shear strength 
than the transverse normal stress. In the cases of [(i0)/(±0)/9O°/9O°] g 
graphite-epoxy, crack surfaces are open and delamination growth is controlled 
by all of the three-dimensional, mixed-mode stress intensity factors, Kj, 

and ^xil* 

(7) While the govern local deformation and delamination fracture, the 
G^ are related to the global structural response and less sensitive to the 
local deformation and fracture. For example, the total G and Gg differ only 
slightly between the cases of a delamination with a closed and an open crack 
tip in the [9/-0/-0/0] graphite-epoxy composite. Thus, may provide a more 
sensitive measure and, thus, better fracture parameters than G^ for evaluating 
the composite delamination growth. 

(8) Influences of lamination variables such as fiber orientation and ply 
thickness on the delamination behavior are significant. Changing fiber 
orientation generally alters failure modes appreciably. For instance, in 


[ (d=0 ) / (±0 )/90°/90°] s graphite-epoxy, the opening-mode (Kj) dominated 
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delamination fracture changes into a shearing-mode (Kjj) controlled failure, 
as the 0 becomes greater than 40°. Also, increasing ply thickness hj/W in 
[ 0/ — 0/ — 0/ 9] graphite-epoxy affects the failure mode from a shearing-dominated, 
open delamination fracture to a closed one. 

(9) Stress intensity solutions and energy release rates appear to be 
independent of crack length during the growth of an edge delamination, as long 
as the crack exceeds a length of few ply thicknesses [3,19]. This unique 
feature in composite edge-delamination is being used for evaluating of 
interlaminar fracture toughness of fiber composite laminates [16,17], provided 
that the and can be calculated accurately by using advanced analytical 
methods such as the present one. 
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TABLE 1 

Delamination Stress Intensity Factors Associated with 
Finite-Element Mesh Refinement in [45°/-45 0 /-45 0 /45°] Graphite-Epoxy 


No. of 
Elements 

Element Size 
(Ah/h) 

* 

K I 

K II 

K III 

72 

0.25 

-0.6458E-1 

0.3424E-2 

-0.7114E 0 

84 

0.125 

-0.6219E-1 

0.2463E-2 

-0.6967E 0 

96 

0.0625 

-0.6141E-1 

0.1796E-2 

-0.6894E 0 

108 

0.03125 

-0 . 6152E-1 

0.1401E-2 

-0.6855E 0 


are normalized by 10 6 e (psi-/in) . 

* 

Closed crack tip with c/a = 0.34. 
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TABLE 2 

Delamination Energy Release Rates G-j^ Associated with 
Finite-Element Mesh Refinement in [45° /-45 0 /-45°/45° ] Graphite-Epoxy 


No . of 
Elements 

Element Size 
(Ah/h) 

* 

G 1 

G 2 

G 3 

72 

0.25 

0 . 

0.6855E-6 

0.1533E 0 

84 

0.125 

0 . 

0.5528E-6 

0.1519E 0 

96 

0.0625 

0 . 

0.4619E-6 

0.1512E 0 

108 

0.03125 

0 . 

0.3943E-6 

0.1508E 0 


^G^ are normalized by 10 6 e 2 (lb-in/in 2 ). 
* 

Closed crack tip with c/a =0.34 
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TABLE 3 

Influence of Fiber Orientation on Crack-Tip Closure Length 
for a Delamination in [0/-0/-0/0] Graphite-Epoxy Composite T 


0 

c/a 

0° 

0.0 

15° 

0.2292 

30° 

0.2834 

45° 

0.3402 

60° 

0.3603 


i 

h- = h ® 0.25 in., b = 2 in.; 
a = 0.25 in. 

& 

at strain e = 10 -6 
o 
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TABLE 4 

Effects of Fiber Orientations on Stress Intensity Solutions 
for Delaraination in [0/-0/-Q/0] Graphite-Epoxy Composite^ 


0 

K i 

K II 

K III 

0° 

0 . 

0 . 

0. 

15° 

-0 . 3782E-1 

0 . 3866E-2 

-0.2379E 1 

30° 

-0.1072E 0 

0 • 8744E-2 

-0.1943E 1 

45° 

-0.6152E 0 

0 . 1401E-2 

-0.6855E 0 

O' 

0 

0 

-0 . 8599E-2 

0.2034E-3 

-0.9270E-1 


+ 

hi = V 12 = 0.25 in., b = 2.0 in.; a = 0.25 in. 
are scaled by lO^e^ [psi-^in] . 



39 


TABLE 5 


Effects of Fiber Orientation on Strain Energy Release Rates 
for a Delamination in [0/-0/-0/0] Graphite-Epoxy Composite 


* 

t 


0 

« 1 tt 

G 2 

G 3 

0° 

0. 

0.0 

0.0 

15° 

0. 

0.6059E-5 

0.2111E 1 

30° 

0. 

0 . 2748E-4 

0.1074E 1 

45° 

0. 

0 . 3943E-6 

0.1508E 0 

o 

0 

0. 

0.1425E-7 

0 . 3597E-2 


i* 

= 0.25 in., b = 2 in.; a = 0.25 in. 

G. are scaled by ( 10 6 e 2 ) [lb-in/in 2 ] . 

1 o 

tt 


is identically zero due to crack-tip closure. 
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9. LIST OF FIGURE CAPTIONS 


Fig. 1 

Fig. 2 
Fig. 3 

Fig. 4 
Fig. 5 
Fig. 6 
Fig. 7 
Fig. 8 
Fig. 9 


Fig. 10 


Fig. 11 


(a) Six-Node Quasi Three-Dimensional Crack-Tip Element with Singular 
Derivatives; (b) Eight-Node Quasi Three-Dimensional (nonsingular) 
Isoparametric Element. 

Arrangement of Special Crack-Tip Elements for the Composite 
Delamination Problem. 
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APPENDIX 1 

Shape Functions N ( p , £ ; 6 ) for the Singular Delamination Crack— Tip Element 

M S 

Corresponding to a singular eigenvalue <5, the shape functions N in local 
triangular polar coordinates (p,£) for the six-node conforming crack-tip 

element in the delamination problem are shown [4] to have the following 
expressions : 

N ii = 1 + K2-2' <S )p-P 6 + 1 ]/(2" 6 -1), (Al-1) 

N i(3+i) = 2" 6 p(l-0/(2" 6 -l) - p 6 + 1 [2S(l-0+(l-0/(2 - 6 -l)], (Al-2) 

N i( 6 +i) = 2~ < 5 pC/(2“ < 5 -1) - p 6 + 1 [2C(1-5)+5/(2" 6 -1)], (Al-3) 

N i(9+i) = " t2p(l-0+2p 6 + 1 (l-C)]/(2" 6 -l), (Al-4) 

N i( 12 +i) = 4 P <S+1 ^ 1 "?)» (Al-5) 

N i(15+1) = “ 2[pe-p 6 + 1 S]/(2" 6 -l), (i » 1,2,3), (Al- 6 ) 

where the singular eigenvalues 6 's are determined in accordance with the local 
crack-surface boundary conditions discussed in Section 2.1, and the p and £ 
are related to the global coordinates by 

(x 9 -x.) tan <p -(y 9 -y.) 

5 = 7 — — — ^ — 7 s-r — ~r* (Ai-7) 

(y 3 -y 2 ) -(x 3 -x 2 ) tan <j>’ 

P = r/f(5), (Al- 8 ) 

in which 

f( 5 ) = {(x 2 -x 1 ) 2 +(y 2 -y 1 ) 2 + 25 [(x 2 ~x 1 )(x 3 -x 2 )+(y 2 -y 1 )(y 3 -y 2 )] 
+ 5 2 [(x 3 ^c 2 ) 2 +(y 3 -y 2 ) 2 ]} 1/2 . 


(Al-9) 
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APPENDIX 2 

Iteration Scheme for Partially Closed Delamination 

For a delamination with crack surfaces in finite-length contact, the 
following iteration scheme is used for determining contact stress and crack 
closure length: 

(i) Assume an initial contact length Cj, and solve Eq. 21 with the side 
conditions Eqs. 23(a) and 23(b) for contact stress and displacement. 

(ii) Check the solution against preset criterion rif. If the compressive 

stress Fj = (^(-c^tt) gives a = | ( F i~ F 0 )/ F o I > and F o is a 

properly selected sealing factor, for example, F Q is set as 10^ 

(psi), in the present calculation.] such that ri^ < rjf and if the 
displacement field is admissible (i.e., no overlapping or 
interpenetration beyond the contact region), we set the crack closure 
length c = cj and terminate the iteration. 

(iii) If the two constraint conditions are not met, a new contact length C 2 > 
c^ is assumed and the procedure of (i) and (ii) is repeated for F 2 and 
u 2 . 

(iv) If F^ < 0 (i = 1,2) and [r^l > |n^| > nj, the next assumed length C 3 
should be C 3 < cj, and repeat (iii) and (iv). 

(v) If either F^ < 0 and |nj| > | n 2 | > or FjF 2 < 0 and | | > the 
next assumed contact length is set as 

c 3 - c 2 + F 2 ( c 2" c l^( F l“ F 2^* (A2-1) 

The iteration from (i)-(iii) continues until sufficient accuracy is 


reached. 
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(vi) As the constraint conditions are satisfied and the difference between 
assumed contact lengths is |c n+1 -c n | < e c , the iteration is terminated 
with c = V2 | c n +l+c n | • 
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